% ***************************************************************************************
% m_calibration_AD.m 
% estimate alpha1,alpha2 and beta1,beta2 and model inversion 
% ***************************************************************************************
clear all; 
close all;
format shortG;

% ==================== SET FUNDAMENTAL PARAMETERS ========================
gravity = 8;             % gravity coefficient (estimation in step2 gravity)   
rho = 0.66;              % discount factor (after 1955) 
rho50 = 0.66;            % --- (for 1950) 
sigma = rho/gravity;     % shape parameter of idiosyncratic shocks (after 1955)  
sigma50 = rho50/gravity; % --- (for 1950)

mu=0.15;                 % expenditure share on floor spaces (consumers) 
gammaL=0.7;              % cost share on labor (production)
gammaH=0.1;              % cost share on floor spaces (production)
eta=4;                   % floor space supply elasticity 

lbalpha=-0.10; ubalpha=0.30;  alphagridint=0.001; 
lbbeta=-0.10;  ubbeta=0.30;   betagridint=0.001; 
alphagrid=(lbalpha:alphagridint:ubalpha);  
betagrid=(lbbeta:betagridint:ubbeta);      
Ngrid=size(alphagrid,2); 
[m,n] = ndgrid(alphagrid', alphagrid');  alphagridmat = [m(:),n(:)]; 
[m,n] = ndgrid(betagrid', betagrid');    betagridmat = [m(:),n(:)]; 

% ========================= LOAD DATA  =========================================  
datafile = '../../data/processed data/quant/Hiroshima_DATA.csv'; 
CITYDATA = readtable(datafile);  

ID = CITYDATA.geocode;   IDunique = unique(ID); 
if size(ID,1) ~= size(IDunique,1) 
    disp(' >>> CHECK DUPLICATED LOCATION DATA <<<');  pause 
end 
clear IDunique 
N=size(ID,1); Area = CITYDATA.area; GEOGRAPHY=[ID Area]; cbddist=CITYDATA.cbddist; 

Rdata36=CITYDATA.pop36; Rdata45=CITYDATA.pop45bomb; 
Rdata51=CITYDATA.pop51; Rdata55=CITYDATA.pop55; 
Rdata60=CITYDATA.pop60; Rdata65=CITYDATA.pop65;     
Rdata70=CITYDATA.pop70; Rdata75=CITYDATA.pop75; 
Rdata50=CITYDATA.pop50; Rdata53=CITYDATA.pop53;
Rdata57 = 0.6*Rdata55+0.4*Rdata60; 
Rdata66 = 0.8*Rdata65+0.2*Rdata70;

Ldata38=CITYDATA.emp38; Ldata45=CITYDATA.emp45bomb; 
Ldata53=CITYDATA.emp53; Ldata57=CITYDATA.emp57; 
Ldata63=CITYDATA.emp63; Ldata66=CITYDATA.emp66;     
Ldata75=CITYDATA.emp75; 

% Read latitude and longitude (for readability of the csv output in GIS)
lat = CITYDATA.lat; lon = CITYDATA.lon;

Ldata50=Ldata53*(sum(Rdata50)/sum(Rdata53)); 
Ldata55=0.5*Ldata53+0.5*Ldata57;  
Ldata60=0.5.*Ldata57+0.5.*Ldata63; 
Ldata65=(1/3)*Ldata63 + (2/3)*Ldata66; 
Ldata70=(5/9).*Ldata66+(4/9).*Ldata75; 

Rdata45=(sum(Ldata45)/sum(Rdata45)).*Rdata45;  
Rdata50=(sum(Ldata50)/sum(Rdata50)).*Rdata50;
Rdata55=(sum(Ldata55)/sum(Rdata55)).*Rdata55;  
Rdata60=(sum(Ldata60)/sum(Rdata60)).*Rdata60;
Rdata65=(sum(Ldata65)/sum(Rdata65)).*Rdata65;  
Rdata70=(sum(Ldata70)/sum(Rdata70)).*Rdata70;
Rdata75=(sum(Ldata75)/sum(Rdata75)).*Rdata75; 

% data for the post-war economy (every 5 years from 1950);
% NOTE: we use data from 1950 to 1975 for estimation
Rdata=[Rdata50 Rdata55 Rdata60 Rdata65 Rdata70 Rdata75];
Ldata=[Ldata50 Ldata55 Ldata60 Ldata65 Ldata70 Ldata75]; 
totalR = sum(Rdata,1); totalL = sum(Ldata,1);  
Popdata=totalR; 
if min(totalR ~= totalL) == 1   % check total pop ; 
    disp('>>> CHECK TOTAL POPULATION AND EMPLOYMENT ! <<< ');   pause
end 

% all data (every 5 years, 1945-75) 
Rdata_all = [Rdata45 Rdata50 Rdata55 Rdata60 Rdata65 Rdata70 Rdata75];
Ldata_all = [Ldata45 Ldata50 Ldata55 Ldata60 Ldata65 Ldata70 Ldata75]; 

% population/employment density; 
NT=size(Rdata,2);   
Rdensdata = Rdata./repmat(Area,1,NT);  
Ldensdata=Ldata./repmat(Area,1,NT);  

% COMPUTE LOWER BOUNDS for Theta; 
Rthetalb = (Rdata(:,2:NT)-Rdata(:,1:NT-1) < 0);  
Rthetalb = (1 - Rdata(:,2:NT)./Rdata(:,1:NT-1)).*Rthetalb;  
Lthetalb = (Ldata(:,2:NT)-Ldata(:,1:NT-1) < 0);  
Lthetalb = (1 - Ldata(:,2:NT)./Ldata(:,1:NT-1)).*Lthetalb; 
Rthetalb = max(Rthetalb); Lthetalb = max(Lthetalb); 
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;

% ================= LOAD FLOORSPACE PRICES ================================
l_q50=CITYDATA.lnQ_50; q50=exp(l_q50*(4/eta)); 
l_q55=CITYDATA.lnQ_55; q55=exp(l_q55*(4/eta)); 
l_q60=CITYDATA.lnQ_60; q60=exp(l_q60*(4/eta)); 
l_q65=CITYDATA.lnQ_65; q65=exp(l_q65*(4/eta)); 
l_q70=CITYDATA.lnQ_70; q70=exp(l_q70*(4/eta));
l_q75=CITYDATA.lnQ_75; q75=exp(l_q75*(4/eta)); 
Qdata=[q50 q55 q60 q65 q70 q75]; % Floor space price matrix  

% ================= LOAD COMMUTING COST ===================================
load('commuting_cost.mat', 'Kappa_mat');  

% ================= SET Theta VECTOR ======================================
thetavec=[0.53; 0.53; 0.53; 0.53; 0.53];

if min(thetavec - thetalb')<0
    disp('>>> CHECK THETA VALUES (LOWER BOUND) ! <<< ');   pause
end 

%  =============  GROUP TRACTS FOR GMM BY DISTANCE FROM CBD ===============
G=5;
[Gmat] = gmmiv(cbddist,G,N);  
GG = sum(Gmat,2); 

%  ========================================================================
%  ==================  MAIN INVERSION OF THE MODEL   ======================
%  ========================================================================
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
Xmat=zeros(N,NT); 
Ymat=zeros(N,NT); 
Kmat = zeros(N,N,NT); 
LQmat=zeros(N,NT); 

Kmat(:,:,NT)=Kappa_mat(:,:,NT).^(-rho/sigma);   % last period 
for t=1:NT-1
    tt=NT-t; 
    theta=thetavec(tt); 
    if t < NT-1; Kmat(:,:,tt)=(Kappa_mat(:,:,tt).^(-rho/sigma)).*(Kmat(:,:,tt+1).^(rho*(1-theta))); 
    else;        Kmat(:,:,tt)=(Kappa_mat(:,:,tt).^(-rho50/sigma50)).*(Kmat(:,:,tt+1).^(sigma*rho50*(1-theta)/sigma50)); end
end 

LQmat(:,NT)=Qdata(:,NT);
for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    if t < NT-1; LQmat(:,tt)=Qdata(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta))); 
    else;        LQmat(:,tt)=Qdata(:,tt).*(LQmat(:,tt+1).^(rho50*(1-theta))); end
end 

K50=Kmat(:,:,1);
Q50=LQmat(:,1); 

%  =================== INVERSION OF FUNDAMENTALS =========================== 
for t=1:NT-1
    tt = NT-t;    % compute from the last period (1975,1970,1965,1960,1955)
    theta=thetavec(tt); 
    K = Kmat(:,:,tt+1); 
    Q=LQmat(:,tt+1);
    Rpre=Rdata(:,tt); Rpost=Rdata(:,tt+1); 
    Lpre=Ldata(:,tt); Lpost=Ldata(:,tt+1);
       
    X0=ones(N,1);  Y0=ones(N,1);     
    Rsyst = @(x)fReval2(x,K,Q,Rpre,Rpost,Lpre,Lpost,rho,sigma,theta,mu,gammaL,gammaH);
    Lsyst = @(y)fLeval2(y,K,Q,Rpre,Rpost,Lpre,Lpost,rho,sigma,theta,mu,gammaL,gammaH);  
    [x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0); X=abs(x_fsolve); X=X./geomean(X);  
    [y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); Y=abs(y_fsolve); Y=Y./geomean(Y);
    Xmat(:,tt+1)=X;    Ymat(:,tt+1)=Y;
end 

X55 = Xmat(:,2); Y55 = Ymat(:,2); 
X60 = Xmat(:,3); Y60 = Ymat(:,3); 
X65 = Xmat(:,4); Y65 = Ymat(:,4); 
X70 = Xmat(:,5); Y70 = Ymat(:,5); 
X75 = Xmat(:,6); Y75 = Ymat(:,6); 

%  ===================== DECOMPOSE FUNDAMENTALS ===========================
amat3D=zeros(N,NT-1,Ngrid^2); bmat3D=zeros(N,NT-1,Ngrid^2);
a_errmat3D = zeros(N,NT-1,Ngrid^2); b_errmat3D=zeros(N,NT-1,Ngrid^2);
for ii=1:Ngrid^2
    alphacomb = alphagridmat(ii,:); alpha1=alphacomb(1); alpha2=alphacomb(2); 
    [a,a_err] = ffund(alpha1,alpha2,Ymat,Ldensdata,rho,sigma,thetavec,N,NT);

    betacomb = betagridmat(ii,:); beta1=betacomb(1); beta2=betacomb(2); 
    [b,b_err] = ffund(beta1,beta2,Xmat,Rdensdata,rho,sigma,thetavec,N,NT) ; 

    amat3D(:,:,ii)=a; a_errmat3D(:,:,ii)=a_err;  
    bmat3D(:,:,ii)=b; b_errmat3D(:,:,ii)=b_err; 
end

%  ===========  GRID SEARCH: 1st Step GMM OBJECTIVE FUNCTION ===============
objvala = zeros(Ngrid^2,1);  objvalb = zeros(Ngrid^2,1);
d_a_err = a_errmat3D(:,2:NT-1,:) - a_errmat3D(:,1:NT-2,:);   
d_b_err = b_errmat3D(:,2:NT-1,:) - b_errmat3D(:,1:NT-2,:);
for ii=1:Ngrid^2
    objvala(ii)=fobjm(Gmat,GG,d_a_err(:,:,ii),N,NT,G,eye(G*(NT-2))); 
    objvalb(ii)=fobjm(Gmat,GG,d_b_err(:,:,ii),N,NT,G,eye(G*(NT-2)));
end
ii_mina = find(objvala==min(objvala)); ii_minb=find(objvalb==min(objvalb)); 

objmatA = reshape(objvala,Ngrid, Ngrid); objmatB = reshape(objvalb,Ngrid, Ngrid); 
[minrowa,mincola]=find(objmatA==min(min(objmatA))); [minrowb,mincolb]=find(objmatB==min(min(objmatB))); 
alpha1est1st=alphagrid(minrowa); alpha2est1st=alphagrid(mincola); 
beta1est1st=betagrid(minrowb);   beta2est1st=betagrid(mincolb);
[alpha1est1st, alpha2est1st, beta1est1st, beta2est1st]

% ======================================================================
% ==========================  2nd Step GMM   ===========================
% ======================================================================
% Change in log residual fundamentals from first step, N x (NT-2) ; 
d_a_err1st = d_a_err(:,:,ii_mina); d_b_err1st=d_b_err(:,:,ii_minb); 

mmadist1st=zeros(G,N,NT-2); mmbdist1st=zeros(G,N,NT-2);
for t=1:NT-2
    mmadist1st(:,:,t)=Gmat.*repmat(d_a_err1st(:,t)',G,1);  mmbdist1st(:,:,t)=Gmat.*repmat(d_b_err1st(:,t)',G,1);   
end
va = zeros(G*(NT-2),G*(NT-2)); vb=va; 
for i=1:N
    mmadistv = reshape(squeeze(mmadist1st(:,i,:)), G*(NT-2),1);  
    mmbdistv = reshape(squeeze(mmbdist1st(:,i,:)), G*(NT-2),1);  
    va = va + mmadistv*mmadistv'./N;   vb = vb + mmbdistv*mmbdistv'./N; 
end

objvala2=zeros(Ngrid^2,1);  objvalb2=zeros(Ngrid^2,1);
wa = inv(va); wb=inv(vb);    % weight matrix 
for ii=1:Ngrid^2
    objvala2(ii)=fobjm(Gmat,GG,d_a_err(:,:,ii),N,NT,G,wa);
    objvalb2(ii)=fobjm(Gmat,GG,d_b_err(:,:,ii),N,NT,G,wb); 
end
objmatA2 = reshape(objvala2,Ngrid, Ngrid); objmatB2 = reshape(objvalb2,Ngrid, Ngrid); 
[minrowa2,mincola2] = find(objmatA2==min(min(objmatA2))); 
[minrowb2,mincolb2] = find(objmatB2==min(min(objmatB2))); 
alpha1est=alphagrid(minrowa2); alpha2est=alphagrid(mincola2); beta1est=betagrid(minrowb2); beta2est=betagrid(mincolb2);
[alpha1est alpha2est beta1est beta2est] 

h_a1 = max(1e-6, sqrt(eps)) * max(1, abs(alpha1est));
h_a2 = max(1e-6, sqrt(eps)) * max(1, abs(alpha2est));
h_b1 = max(1e-6, sqrt(eps)) * max(1, abs(beta1est));
h_b2 = max(1e-6, sqrt(eps)) * max(1, abs(beta2est));

alpha1est_se = alpha1est + h_a1;  
alpha2est_se = alpha2est + h_a2;
beta1est_se = beta1est + h_b1;    
beta2est_se = beta2est + h_b2; 

% compute moment function for parameter perturbation: 
[~,a_err_se1] = ffund(alpha1est_se, alpha2est,Ymat,Ldensdata,rho,sigma,thetavec,N,NT) ;  % N x (NT-1): period 2,3,...,T
[~,a_err_se2] = ffund(alpha1est, alpha2est_se,Ymat,Ldensdata,rho,sigma,thetavec,N,NT) ;  
[~,b_err_se1] = ffund(beta1est_se, beta2est,Xmat,Rdensdata,rho,sigma,thetavec,N,NT) ;    
[~,b_err_se2] = ffund(beta1est, beta2est_se,Xmat,Rdensdata,rho,sigma,thetavec,N,NT) ;    
d_a_err_se1=a_err_se1(:,2:NT-1)-a_err_se1(:,1:NT-2);   d_a_err_se2=a_err_se2(:,2:NT-1)-a_err_se2(:,1:NT-2); 
d_b_err_se1=b_err_se1(:,2:NT-1)-b_err_se1(:,1:NT-2);   d_b_err_se2=b_err_se2(:,2:NT-1)-b_err_se2(:,1:NT-2);  

mma_se=zeros(G,N,NT-2,2); mmb_se=zeros(G,N,NT-2,2); % G x N x NT-2 x 2 
for j=1:2
for t=1:NT-2
    if j==1
        mmat_se=Gmat.*repmat(d_a_err_se1(:,t)',G,1);   mmbt_se=Gmat.*repmat(d_b_err_se1(:,t)',G,1);   % G x N 
    else    
        mmat_se=Gmat.*repmat(d_a_err_se2(:,t)',G,1);   mmbt_se=Gmat.*repmat(d_b_err_se2(:,t)',G,1);  
    end 
    mma_se(:,:,t,j) = mmat_se;   mmb_se(:,:,t,j) = mmbt_se;
end
end 

% disp(' >>> COMPUTE MOMENT FUNCTION FOR ESTIMATED PARAMETERS <<< ')
[aa,a_err] = ffund(alpha1est,alpha2est,Ymat,Ldensdata,rho,sigma,thetavec,N,NT) ;  
[bb,b_err] = ffund(beta1est,beta2est,Xmat,Rdensdata,rho,sigma,thetavec,N,NT) ;   
d_a_err=a_err(:,2:NT-1)-a_err(:,1:NT-2);    d_b_err=b_err(:,2:NT-1)-b_err(:,1:NT-2);    
mma=zeros(G,N,NT-2); mmb=zeros(G,N,NT-2);% G x N x NT-2
for t=1:NT-2
    mma(:,:,t)=Gmat.*repmat(d_a_err(:,t)',G,1); 
    mmb(:,:,t)=Gmat.*repmat(d_b_err(:,t)',G,1);   % G x N
end

grad_a1 = squeeze(sum(squeeze(mma_se(:,:,:,1)) - mma, 2));   
grad_a2 = squeeze(sum(squeeze(mma_se(:,:,:,2)) - mma, 2));   

grad_b1 = squeeze(sum(squeeze(mmb_se(:,:,:,1)) - mmb, 2));   
grad_b2 = squeeze(sum(squeeze(mmb_se(:,:,:,2)) - mmb, 2));   

grad_a1 = reshape(grad_a1, G*(NT-2), 1) / h_a1 / N;   % M x 1, M = G*(NT-2)
grad_a2 = reshape(grad_a2, G*(NT-2), 1) / h_a2 / N;
grad_a = [grad_a1, grad_a2]; 

grad_b1 = reshape(grad_b1, G*(NT-2), 1) / h_b1 / N;
grad_b2 = reshape(grad_b2, G*(NT-2), 1) / h_b2 / N;
grad_b = [grad_b1, grad_b2]; 

va = (va + va.')/2;
vb = (vb + vb.')/2;

s_a = grad_a' * (va \ grad_a);     
s_b = grad_b' * (vb \ grad_b);     
v_alpha_gmm = sqrt(1 ./ (s_a * N));
v_beta_gmm  = sqrt(1 ./ (s_b * N));

% =================== SAVE PARAMETER RESULTS ========================
objresult=[gravity,rho,rho50,sigma,sigma50,alpha1est,alpha2est,beta1est,beta2est,v_alpha_gmm(1,1),v_alpha_gmm(2,2),v_beta_gmm(1,1),v_beta_gmm(2,2)];
resultparam=array2table(objresult,'VariableNames',{'gravity','rho','rho50','sigma', 'sigma50', ... 
    'alpha1','alpha2','beta1','beta2','sdalpha1','sdalpha2','sdbeta1','sdbeta2'}); 
writetable(resultparam,'../../output/quant/table_parameter_estimates_AD.csv'); 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                       FUNCTIONS
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------- FUNCTION FOR CONSTRUCTING IV --------- 
function[gmat]=gmmiv(ivvar,nG,N); 
% ivvar: variable for constructing IV
Ground=round(N/nG); gn=zeros(N,1); gmat=zeros(nG,N);
[~,sortn] = sort(ivvar,'ascend'); 
Gthres=[Ground+1; 2*Ground+1; 3*Ground+1; 4*Ground+1; 5*Ground+1; 6*Ground+1; 7*Ground+1; 8*Ground+1; 9*Ground+1];
for i=1:N
    ii=find(sortn==i);
    if ii<Gthres(1); gn(i)=1; 
    elseif ii >= Gthres(nG-1); gn(i)=nG; 
    else; 
        for g=2:nG-1 
        if Gthres(g-1) <= ii && ii < Gthres(g); gn(i)=g;  
        else; continue  
        end 
        end 
    end 
end
for i=1:N
    gi = gn(i);
    gmat(gi,i)=1;
end
end 

% --------  FUNCTION FOR COMPUTING CHARACTERISTICS BASED ON GRAVITY STURCTURE --------
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)   % eq D.16. X corresponds to Theta

N = size(DRpre,1); 
X=abs(X);  % avoid negative values
X=X./geomean(X);
%xtemp= K.*(repmat(X',N,1).^(rho/sigma)); 
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 

Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 

Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) % eq D.17. Y corersponds to Omega

N = size(DRpre,1); 
Y=abs(Y);  % avoid negative values
Y=Y./geomean(Y);
%ytemp= K.*(repmat(Y,1,N).^(rho/sigma)); 
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 

Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 

Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end

% --------  FUNCTION FOR COMPUTING THE EXOGENOUS FUNDAMENTALS --------
function[c, c_err] = ffund(param1,param2,Zmat,densdata,rho,sigma,thetavec,N,NT)  

Z0=Zmat(:,NT); cmat=zeros(N,NT-1);
for t=1:NT-1
    denspost=densdata(:,NT-t+1); denspre=densdata(:,NT-t); 
    invspill = (denspost.^(-param1)).*(denspre.^(-param2));
    cc = Z0.*invspill;
    cmat(:,NT-t)= cc; 
    theta=thetavec(NT-t); 
    Z0=Zmat(:,NT-t).*(Z0.^(-rho*(1-theta)));     
end

c = log(cmat);   % N x NT-1
c_te = sum(c,1)./N;
c_err = c-repmat(c_te,N,1); 
end 

% ----------  FUNCTION FOR OBJECTIVE FUNCTION ------------
% Gmat: G (number of grids) x N (number of locations)  
% cmat: N (number of locations) x NT-2 (number of periods well defined for changes)
% Wmat: weight matrix, G*NT-2 (number of grids) x G*NT-2 
% objdist: scalar 
function[objdist] = fobjm(Gmat,GG,cmat,N,NT,G,Wmat); 
mdist = zeros(G,N,NT-2);
for t=1:NT-2
    mmt = Gmat.*repmat(cmat(:,t)',G,1);  % G x N
    mdist(:,:,t) = mmt;                 
end
mdist = squeeze(sum(mdist,2));   %  G x NT-2
mdist = reshape(mdist, G*(NT-2),1); 
objdist=(mdist./N)'*Wmat*(mdist./N);
end
